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Abstract. We develop a generic method to compute the dynamics induced by 
quenches in completely connected quantum systems. These models are expected 
to provide a mean-field description at least of the short time dynamics of finite 
dimensional system. We apply our method to the Bose-Hubbard model, to a 
generalized Jaynes-Cummings model, and to the Ising model in a transverse field. We 
find that the quantum evolution can be mapped onto a classical effective dynamics, 
which involves only a few intensive obscrvables. For some special parameters of the 
quench, peculiar dynamical transitions occur. They result from singularities of the 
classical effective dynamics and are reminiscent of the transition recently found in the 
fermionic Hubbard model. Finally, we discuss the generality of our results and possible 
extensions. 
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1. Introduction 

Thanks to the fast experimental progress on cold atoms [H EJ [3], where direct control 
over the parameters of the effective Hamiltonian is available, many theoretical questions 
about out of equilibrium dynamics of quantum systems have been raised and started 
to be addressed. Among them, we cite the problem of thermalization for an isolated 
system (see e.g. [U Ej) the existence of long-lived out of equilibrium states (see e.g. 
[EJ [3 El E]), dynamical transitions out of equilibrium [TUJ EJ [121 HB], etc. 

In this work we focus on the off-equilibrium dynamics induced by quantum 
quenches, i.e. the evolution of an isolated quantum system after a sudden change 
of a parameter of the Hamiltonian. This problem has been intensively studied these 
last years. It provides an useful idealization of the protocols followed in experiments, 
where the change of parameters takes instead place at finite rates. The literature on 
quantum quenches is already quite broad, see for example the reviews [T^| [T5| IT6] . One 
dimensional systems have been intensively studied by exact analytical and numerical 
methods. Results for higher dimensional systems are instead scarcer; the previous 
methods cannot be applied and one has to resort to approximations of some kind. The 
ones that have been more used are the Bogoliubov method [17] . path integral saddle- 
point expansions [T8| [19], the large- N limit [20] and mean-field (or large dimension) 
approximations [21} \TU\ HH CEJ [22] . Interestingly, mean- field approaches revealed that 
quantum quenches may lead to out of equilibrium dynamical transitions. This was first 
found in the analysis of the Fermionic Hubbard model by time dependent Dynamical 
Mean Field Theory (t-DMFT) in [10]. Later, a confirmation and explanation of this 
result was found by using a Gutzwiller time-dependent Ansatz in [IT]. In [12] we 
solved the Bose Hubbard model on a completely connected lattice and found a very 
similar dynamical transition at integer fillings. This was also noticed in the hard-core 
Bose-Hubbard model in a superlattice potential [23]. Finally, clues to argue that these 
dynamical transitions are generic even beyond mean-field were recently presented in [13] 
by mapping the problem to the one of classical phase transitions in films. 

In this article we first describe a generic approach to solve the dynamics of 
quenches in completely connected quantum models. Then, we apply our method 
to the Bose-Hubbard model, the generalized Jaynes-Cummings model and the Ising 
model in a transverse field. We find that the dynamical transition, discovered for 
the Hubbard model, is a systematic dynamical effect present in completely connected 
systems characterized by a quantum phase transition in equilibrium. A shorter version 
of this work, that focused only on the Bose-Hubbard model, appeared in [T2] . 

2. Summary of the method and results 

In this article, we focus on the out of equilibrium dynamics in several completely 
connected models: the Bose Hubbard model (BHM), the transverse field Ising model 
(IM) and the Jaynes-Cumming model (JCM). The reason for focusing on systems 
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defined on completely connected graphs is that this allows one to study a well defined 
model and — at the same time — to obtain an approximate solution for finite dimensional 
lattices. Indeed, such models are related to several approximations used in the literature. 
Among the most known ones, we cite the limit of infinite dimensions for models defined 
on a hyper-cubic d- dimensional lattice^ and the Gutzwiller Ansatz, a widespread method 
to obtain mean field equations. Actually, in the case of completely connected models, 
the Gutzwiller Ansatz can be shown to be exact by using a Hubbard-Stratonovich 
transformation. 

We shall focus on out of equilibrium dynamics induced by a quantum quench. This 
procedure is defined as follows. Let H(\) be the quantum Hamiltonian, which depends 
on a coupling A. The system is prepared in its ground state \ip(t < 0)) = |GS(Aj)) at 
some coupling Aj. At t > 0, the system is driven out of equilibrium in a controlled 
fashion, tuning the coupling X(t) in time according to a predefined procedure. Except 
for quasistatic procedures, often called "adiabatic" , the wave function becomes different 
from the ground state \tj)(t > 0)) ^ |GS(A(t))). We call sudden quench the procedure in 
which the coupling is suddenly switched to a final value A/: X(t) = (Xf — Aj) 6{t) + Aj. 

As we shall show, the quantum dynamics in a completely connected system can be 
solved by mapping the unitary evolution onto an effective model undergoing Newtonian 
dynamics. This is a drastic simplification which is possible thanks to the symmetry 
of the completely connected Hamiltonian under any permutation of sites. We restrict 
our analysis to cases where \ip{t = 0)) is the ground state at some coupling Aj. Thus, 
\i/j(t = 0)) is also symmetric under permutation of sites and since both the initial state 
and the Hamiltonian are symmetric, the generic unitary evolution takes place in the 
sector of symmetric states only. In this subspace, the states are parametrized using 
a few local macroscopic observables. The unitary evolution can then be written as a 
Schrodinger equation in the symmetric space, which involve an effective H = V' 1 , where 
V is the number of sites of the system. Thanks to this property, at the thermodynamic 
limit, the entire dynamics of the system can be encoded in the one of few macroscopic 
variables, which undergo an effective classical Hamiltonian evolution. 

We first apply this approach to the Bose-Hubbard model, our main interest because 
of its applicability to cold atom experiments. In this case, it is the onsite repulsion U 
that plays the role of the coupling A. In a first stage we make the additional assumption 
that the number of bosons per site is less or equal to n™ ax = 2. The average superfluid 
order after the quench is a non-monotonous function of \Uf — Ui\. Actually, it decays 
logarithmically to zero at some special values of [/, and Uf. For these special quenches, 
the superfluid order relaxes exponentially to zero. We call this peculiar feature a 
dynamical transition. Surprisingly, the microcanonical equilibrium characterized by 
same energy, towards which the system would relax on large times if it were able to 



thermalizq§l, has non-zero superfluid order. Thus, this transition is a purely dynamical 
J This consists in a technique similar to the ones developed in [24] . 

§ Thcrmalization is very likely to occur in finite dimensions for the Bose-Hubbard model, because the 
system is not integrable. Instead for completely connected models we do not expect thcrmalization 
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phenomenon, which has nothing to do with the fact that the superfluid order vanishes 
in the high temperature phase. We extend our analysis to n™ ax > 3. In this case the 
effective dynamics involves n™ ax — 1 degrees of freedom, and the classical trajectories 
can be either regular or chaotic. Apart from this difference, the dynamical transition is 
found to be qualitatively unchanged and the quantitative differences with n™ ax = 2 are 
small. 

Finally, we focus on two other systems, for which some out of equilibrium properties 
have already been studied, a generalized Jaynes-Cummings model |25J, and the Ising 
model in a transverse field [21]. We apply our formalism to both systems and show that 
also in these cases there is a dynamical transition for sudden quenches. 

This paper is organized as follows: In section (3[ we briefly describe the three 
considered models. In Section HI we describe how the effective classical dynamics can be 
derived for an arbitrary Hamiltonian. Then in section [5j we derive the classical dynamics 
for the Bose-Hubbard model with truncation n™ ax = 2, and analyze quenches and the 
dynamical transition in section The general case n™ ax > 3 is considered in section 
[7J The effective dynamics and quench properties of the generalized Jaynes-Cummings 
model are described in section [HI and in section for the Ising model. Section [TU] and 
ITTI contain discussions and possible extensions of our work. Appendix A is devoted 
to prove that the same effective dynamics is recovered with a Gutzwiller Ansatz wave 
function, and in Appendix B the WKB eigenstates of the completely connected model 
are discussed. 



3. Definition of the models 



3.1. Bose-Hubbard model and the Mott- Superfluid transition 

The Bose-Hubbard lattice model has regained a lot of interest recently since it can 
be realized and studied in experiments on cold atoms [26J. The Hamiltonian of its 
completely connected version, suited to describe the limit d — > oo of the lattice model, 
reads 

# = ^£ni(ni-l)-^£&}k (1) 

where b\, bi are the bosonic creation and annihilation operators, satisfying [£>*,&,•] = 
and rij = b\bi. The first term is an on-site repulsion between bosons, the second is a 
tunneling term, of amplitude J, and rescaled by the volume V (number of sites). This 
rescaling is needed to obtain a well-defined thermodynamic limit V —t oo. The limit of 
infinite dimensions that is equivalent to this model is obtained by taking a coupling ^ 
for a BHM defined on a d- dimensional hyper-cubic lattice. 

This model undergoes a quantum phase transition at commensurate fillings: (hi) = 
n with n integer. See [27] for a detailed discussion of the phase diagram. Notice that 

to the Gibbs ensemble, as we shall discuss later. This implies that in the large dimension limit the 
timescale for equilibration diverges as a function of d. 
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since we consider the dynamical behavior, the number of particles is fixed because it 
corresponds to a quantity conserved by the dynamics. In consequence we shall discuss 
the phase diagram in terms of the density and not the chemical potential, as it is instead 
done usually. Let us focus first on cases where the number of particles per site is an 
integer. To grasp why there is a quantum phase transition, it is helpful to compare the 
ground state in the two limits U — > oo (or J = 0) and U — > 0. In the first case, the 
ground state is diagonal in the occupation number = (8>i|rij = n). It remains the 
same for all U < U c , and is called a Mott insulator. It is incompressible because adding 
or removing a boson requires an energy of the order of U. In the second case, U — 0, the 
ground state is akin tqjjj a product of coherent states \ip) = C Yli e ~ abi |0)- For U — 0, 
and also moderate values of U, the ground state is superfluid and compressible. The 
usual order parameter {bi) is zero because the density is fixed, hence the appropriate 
parameter is the off-diagonal long range order measured by |^o| 2 = lim c z ij - 5 .oo(^i) • In 
the completely connected model, all distances dij between two different sites are equal 
to 1. This is the largest distance in the problem. In consequence one can defin^jl: 

|* | 2 = (^> i^j (2) 

Because the system is completely connected, phonons are absent and the spectrum 
always has a gap. Beyond a critical coupling U c , the order parameter \^o\ 2 vanishes 
and the ground state becomes a Mott insulator. This is the Mott-superfluid quantum 
phase transition. For non-integer fillings, the system is instead always superfluid and 
compressible as it can be understood considering perturbation around the U — > oo limit. 
At non zero temperature, there is a second order phase transition from the superfluid 
phase to a Bose gas at a critical temperature T C (U). 

3.2. Generalized Jaynes-Cummings model and superradiance transition 

The second model that we consider is the generalized Jaynes-Cummings model studied 
in [25], which describes N distinguishable two- level systems in interaction with a 
single quantized electromagnetic mode. It is useful in various contexts, and has 
been suggested as a qualitative description of the BEC/BCS crossover in fermionic 
condensates, of molecular magnetism and of the formation of dimers by pairing of two 
bosons. Recently, new experiments in the setting of cavity quantum electrodynamics 
[28] are realizations of the Dicke Hamiltonian [29], which can be mapped onto the 
Jaynes-Cummings Hamiltonian if the "rotating wave approximation" is made. The 
dynamics of sweeps, starting from an empty bosonic mode, was studied in [25] using 
quasiclassical approximations and the truncated Wigner approximation. Here, we 
describe the dynamics of sudden quenches from the broken symmetry phase. 

|| Formally, this state is the true ground state only in the grand canonical ensemble. However, the 
canonical and grand canonical ensemble are equivalent up to 1/V corrections. 

% In equilibrium, it is easy to check that the definition below gives back the usual one used in the 
grand-canonical ensemble, \^q\ 2 = \(b)\ 2 
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The model is defined as follows. There is one bosonic [b, b'] = 1 degree of freedom — 
the electromagnetic mode — and one SU(2) spin S = N/2. The spin degree of freedom 
keeps track of the number of excited two-level systems, n e xc = S z + N/2. With n = ¥b, 
the Hamiltonian is 

H = u S z + un + -j= (b ] S~ + bS + ) (3) 



It includes a potential energy for bosons (cu), for the spin (u ), and a coupling g between 
the two. Each time a boson is lost, a two-level system is excited, and reciprocally: the 
quantity Q = S z + tfb is conserved by the dynamics. It is thus possible to define 
2A = Uq — oj such that up to a constant term one finds: 

H = -2\h + ^=(PS~ + bS + ) (4) 



In the following we take g = 1 and consider A in units of g, and the time in units of 
h/g. We consider the regime Q > N/2 (which can be reduced to Q = N/2), for which 
there is a quantum phase transition. At T = 0, the system is in the normal ground 
state (n) = if A < A c = — 1 and is in the super-radiant ground state for A > A c where 
(h) ^ 0. This is the super-radiance quantum phase transition. This transition is the 
consequence of the competition between the two terms in the Hamiltonian: the first 
one favors removing bosons (if A < 0), whereas the second plays the role of a kinetic 
energy and it can lower the energy if the bosons modes are filled, see section [HJ As for 
the Bose-Hubbard model, there is a finite temperature phase transition at T C (A), above 
which the super-radiant phase disappears. Notice that this phase transition is usually 
studied in the canonical/mi crocanonical ensemble, whereas here it is more natural to 
focus on given value of Q since this quantity conserved by the dynamics. 

3. 3. Ising model in a transverse field and the ferromagnetic transition 

The Ising model in a transverse field is a paradigm of quantum phase transitions [30] . 
The Hamiltonian of its completely connected version reads 

£ = -2^£^- r £# (5) 

ij i 

where J is the ferromagnetic coupling and V the transverse field. For simplicity, we 
set J = 1, measure T in units of J and time in units of H/J. The Hamiltonian can be 

written using a single spin S = J2i &i an d reads 

H = — -(S Z ) 2 -TS X (6) 
2N K ' v ; 

The large N limit corresponds to the large spin S = N/2 limit, which is also the classical 

limit. The ground state is found, minimizing the corresponding classical HamiltonianE] 

S — > S = S{sin8 coscj), sin 6 sin 0, cos#} [21]. The minimum is at = 0, and for T < 1/2, 
sin# = 2r, which is a ferromagnetic ground state, in which the spins align toward 

+ The kinetic term of the classical limit is nontrivial to obtain, but we do not need it for the moment. 
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the z-axis. For r > 1/2, = it and the ground state is called quantum paramagnet, 
because it is unoriented in the z-basis |±) i: \tp) = —pv ®% (|+)t + |— }•;)• The quantum 
phase transition from z-ferromagnet to quantum paramagnet takes place at T c = 1/2. 
Starting from the ferromagnetic ground state and increasing the temperature the system 
undergoes a second-order phase transition at a finite temperature T C (T). 

4. Quantum quenches in completely connected models: mapping to an 
effective classical dynamics 

4-1. Site permutation symmetry and classical Hamiltonian dynamics 

In the following, we show how the dynamics induced by quantum quenches, in arbitrary 
completely connected models, can be mapped onto an effective classical Hamiltonian 
dynamics. As already stated, any model defined on a completely connected lattice has 
a Hamiltonian which is symmetric under any permutation of sites. Thus, one expects 
that the ground state of the model is also site permutation symmetricE]. We shall study 
quenches starting from the ground state of the system. Since the site permutation 
symmetry is conserved by the unitary evolution, the dynamics only takes place into 
the subspace of symmetric wave functions. This is a drastic simplification, because 
symmetric states can be described using a few variables only. 

To clarify this point in a concrete but simple case, let us see what these variables 
are for the Bose-Hubbard model, with the additional constraint that there are n™ ax = 2 
or less bosons per site. First, remark that given a particular Fock state the 
linear combination of every possible site permuted Fock state is a site permutation 
symmetric state \S({ni})) ~ J2p \{ n P{i)})- This new state is completely characterized 
by the fraction xq, x±, X2 of sites with 0, 1 and 2 bosons respectively, and we call it 
\xq,Xi,X2) = \x) where x is a shorthand for all variables. Since the Fock states form a 
basis of generic states, the {\x)} states form a basis of the symmetric sector. In order 
to express the Schrodinger evolution in this basis, we have to compute the transitions 
elements (x'\H\x). Typically, only few transitions W m (x) = —h{x + m/V\H\x) are 
allowed^, with m = {mo, mi, . . .} a vector of integers. For example in the Bose-Hubbard 
model, the matrix elements of the Hamiltonian connect states that differs by one boson 
jump. Thus Xq, X\ and x<i can only differ by 1/V or 2/V. For example, after a jump 
|2f ... Oj ...)-> |li .. . lj . . .), x[ = xi + 2/V, x' = xo-l/V and x' 2 = x 2 - 1/V. This 
transition is labeled m = {mo = —l,mi = 2,m>2 = —1}- Moreover, the reverse move is 
a transition characterized by —m, with same amplitude W- m (x) = W m (x) at dominant 
order in V. 

Thanks to the "locality" of transition elements between symmetric states, the 
Schrodinger equation in the symmetric basis = ip x (t)\x) takes a simple form. 

* This symmetry may be spontaneously broken in systems with attractive interactions, which we do 
not consider here. 

(t The minus sign in this definition is for later convenience. 
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The Schrodinger equation projected on (x\ reads (x\id t \ip) = (x\H\ip). Denoting the 
diagonal transition element D(x) = y(x\H\x), we get 



x+m 

m 

= v(D{x)-2^W m {x) co S h{ mi dxjV))^ 



(7) 



where we made use of ip x + m /v = Gxp{Tn,id Xi /V)ij) x with an implicit summation over the 
index i. Strikingly, the Schrodinger equation (J7J) involves an effective H = 1/V, thus 
the regime of interest V — > oo corresponds to the classical regime. Furthermore, for the 
ground state, the initial wave function is a narrow wave packet of width 1/yV (this 
general property is discussed in more detail for the Bose-Hubbard model in section l5\2l) . 
Thus, by the Heisenberg uncertainty principle, we expect that momentum fluctuations 
are proportional to \/h oc y/l/V. Therefore the evolution of the wave-packet can be 
fully described in the thermodynamic limit by its average (or center) x(t) = (x) and 
average momentum p(t) = (p). Both quantities evolve following a classical Hamiltonian 
dynamics obtained from the quantum one by replacing p = ^d x — > p(t) and x — > x(t): 

yd t ^x = (D x -2W x cosh{2d x /V))i> x 

= f D x - 2 Y, m W m(x) COS^iP^Jlpx = H^ x 

H[x,p] = D(x) -2j2W m (x)cos(miPi) (9) 

m 

The classical Hamiltonian evolution of the variables ii(t) = dH/dpi and Pi(t) = 
—dH/ dxi yields the evolution of the wave packet after the quantum quench, and give 
access to all observables as a function of time. A more careful analysis is performed in 
section 15.31 and fully supports this picture. 

As a conclusion, the analysis of quench dynamics in connected models is tractable, 
and takes the form of an effective classical Hamiltonian evolution. The initial condition 
is provided by the ground state values of Xi,pt at the initial coupling f/j. The sudden 
quantum quench dynamics is then described by the classical dynamics of x(t),p(t) 
induced by the effective Hamiltonian H(Uf), with initial conditions Xi,pi. 



5. Effective classical Hamiltonian equations in the Bose-Hubbard model 
with truncation n™ ax = 2 

In the following, we study the Bose-Hubbard model with two or less bosons per site. 
We derive explicitly the effective classical Hamiltonian in section 15.11 The nature of 
the ground state, in particular the peaking of the wave packet when V — > oo and the 
quantum phase transition, are discussed in section 15.21 A more thorough derivation of 
the effective Hamiltonian dynamics is presented in section 15.31 Finally, we conclude by 
describing the phase space of effective trajectories in section l5~4l 
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5.1. Schrddinger equation in the site permutation symmetric basis and effective 
Hamiltonian 

In this section, the method sketched above is applied to the Bose-Hubbard model with 
truncation < n™ ax = 2. This restriction is well-suited to derive analytical expressions, 
and preserves the main features of the model, such as the Mott insulator to superfluid 
quantum phase transition. The case of a larger number of bosons per site is considered 
in section [7J The results remain qualitatively the same. 

The symmetric states are written under the form \xq, £i, x 2 ). Since the total number 
of sites V = V(x +Xi+x 2 ) is fixed, and the overall density of bosons n = N/V = X\+2x 2 
is conserved by the dynamics, the symmetric states can be labeled by one variable only. 
Choosing x\ as such a variable, the symmetric normalized wave functions are denoted 
\x\) = (A r !A r i!A r 2 !/V A !) 1 / 2 Y2{m} I n i5 n 2, ■ ■■,n>v), where Ni = Vx{ and Yl' means the sum 
over Fock states of fixed fraction x\. In the following for simplicity of notation we drop 
the subindex one and use x instead of X\. Note that contrary to the previous paragraph, 
x is now just a number and not a vector. 

To compute the transition rates, we proceed as follows. The repulsion term 
1/2 J2i n i{ n i ~ 1) is diagonal in the \x) basis. The tunneling term —1/V b]bi allows 
transitions from \N , N 1} N 2 ) to three states, \N -1, Nt+2, N 2 -l), \N +1, N x -2, N 2 +l) 
and \Nq, Ni, N 2 ). Using the previous notations, the only possible transitions correspond 
to m = 2, m = —2, m = representing respectively (x + 2/V\H\x), (x — 2/V\H\x) and 
(x\H\x). The first amplitude reads 

(x + 2/V\H\x) = 

/ (AT, - i m H- 2 M N 2 - 1)1 N^m y> ± W }\H±\ {ni] ) < 10 » 

There are V\/(N \Ni\N 2 \) factors on the ket side. For each of these factors, 
there are N N 2 transition to a state with N[ — N± + 2. Therefore using that 
b\bj\0i)\2j) = v^|li)|lj), and after a partial cancellation of the normalization factors, 
we obtain: 

(x + 2/V\H\x) = rpjVoiV 2 (Aq + 1)(N X + 2)] 1 / 2 (11) 

In the following, only dominant contributions of order V are kept, the sub-leading ones 
are dropped since they do not matter for dynamics taking place on times not diverging 
with the system size. Rewriting xq and x 2 as functions of x, the transition rates are 

(x-2/V\H\x) = -Vx[(2-x-n)(n-x)/2}^ 2 =-VW x 
(x\H\x) = VU(n-x)/2-Vx(2 + n-3x)/2 =VD X (12) 

(x + 2/V\H\x) = (x-2/V\H\x) 

The Schrddinger evolution on ip x {t) in the site permutation symmetrical basis is 
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(x\id t \ip) = {x\H\%fj) 

id t ip x = V D x ijj x -VW X (i> x+2 /v + i>x-2/vj 

= v(f) x - 2W x coslx{2d x /V))^ x (13) 
= (^-2^008(2^))^ 

Therefore we have obtained that the Hamiltonian in the subspace of symmetric states 
reads 

H = D X - 2W X cos(2j3) (14) 

Note that the dimension of the Hilbert space has been reduced from e v to V in this 
case. The case n™ ax = 2 is especially convenient because the effective classical motion 
takes place in one dimension, therefore the equation of motion are integrable and the 
evolution easy to understand and describe. 



5.2. Nature of the ground state and quantum phase transition 

As in section HI the Schrodinger equation (TIB"]) involves an effective h = 1/V, and its 
classical limit is obtained through x — > x(t) and p — > p(t). 

H[x,p] = D{x) - 2W(x) cos(2p) (15) 

The possibility of describing the quantum dynamics after a quench in terms of 
classical dynamics relies on the fact that the initial wave function is initially of small 
width (~ 1/y/V). In order to check this, let's first expand the eigenvalue equation at 
lowest order in 1/V: 

Eip x = (D x - 2W X - ^d 2 x )^ x 

Since the kinetic term has a factor 1/V 2 in front, the ground state corresponds to the 
minimum of D x — 2W X and small quantum fluctuations around it. Indeed, around the 
minimum, x GS , a quadratic expansion of W x and D x maps this problem onto a quantum 
harmonic oscillator with m = Wa; GS /8 and uj = \J d x (D — 2W)\ XGS /m. Thus, we find 
that the ground state wavefunction is centered around the absolute minimum x GS of the 
potential D(x) — 2W(x) and has a width o ~ y^h/moo, which is of the order of 1/y/V. 
An appropriate expression for any eigenstate is provided by the WKB expansion is given 



in section Appendix B 



The quantum phase transition taking place as a function of U can be recovered 
within this formalism as we now show. We recall that for the Bose-Hubbard model, 
the quantum phase transition and the Mott phase are present at commensurate fillings 
only, so we have to focus on n = 1 for n™ ax = 2. The ground state of the quantum 
Hamiltonian corresponds to the global minimum of the effective Hamiltonian ( TIB"]) , given 
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by = ^ = 0. These equations lead to p = and x equal to the value x GS , which 
verifies 9 ( D ( X )- 2W ( X )) — g. It is easy to check that 



1 if U > U c , Mott insulator ground state 

~ — < ^ $ U < U c , Superfluid ground state ^ 

where U c — 3 + 2\/2 ~ 5.82843 (note that x < 1). In the thermodynamic limit, the 
Mott insulator ground state is \tiq — 1, ri\ — 1 . . . ) = \x — 1). Instead the superfluid 
ground state is \x — x ), with 1/2 < x < 1. Therefore, in the superfluid state, 
a fraction of sites x-z — (1 — x GS )/2 > contain two bosons. |^o| 2 = ^72 ^l^i) 
is proportional to the intensive kinetic energy, and can be computed using the identity 
|\I/o| 2 = — 1/ J(E — U/2(J2i n i( n i — 1))) = ^gs(1 — x GS )U c /2. It is the order parameter for 
the transition: it is positive in the superfluid phase and vanishes in the Mott insulating 
one. 

5.3. Effective classical evolution of wave packets 

In the following we show in detail how the classical Hamiltonian dynamics emerges from 
the quantum evolution of symmetric states in the thermodynamic limit. For simplicity 
we consider the situation of section [51 where n™ ax < 2, but the argument is more general. 
The Schrodinger equation reads 

^%d t i) x = (p x - 2W X cosh(2<V^))^ (17) 

A ground state wave function is a wave packet characterized by a small width of the 
order of 1/yV as shown in section 15721 The width broadens after a long time, which 
diverges with V. On times that do not diverge in the thermodynamic limit, the ground 
state wave function and its subsequent evolution can be written as 

i/)(x,t) = e - yf ^ t) (18) 

with f(x,t) = g(x,t) — i9(x,t), g and 9 being respectively the envelope and the phase 
of the wave-packet. The envelope g{x,t) has a maximum at x(t), and expanding g(x,t) 
around x(t) shows that indeed this wave function describes a packet of width l/\/V. 

After these preliminary considerations, we proceed and evaluate the evolution of 
the wave function, plugging the expression (1T81) of ip(x,t) into ( |T71) : 

- id t f(x, t) = D x - 2W X cosh(2aj) (19) 

For the amplitude and the phase, this leads to: 

d t g =-2W x sm(2d x 6)smh(2d x g) 

d t e =-D x + 2W x cos{2d x 6)cosh(2d x g) [ ' 

The position of the peak x(t) corresponds to the maximum of the amplitude, thus it 
satisfies the implicit equation d x g\ x ( t ) = 0. To obtain its evolution, we differentiate it 
with respect to time: 

^dlg + d t d x g = (21) 
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The term d t d x g above can be computed from (|2"U|) . At the point where d x g\ x ^ = 0, (1211) 
becomes: 

^ = W x sm(2d x 6)\ x{t) (22) 

This equation involves only one unknown quantity, d x 8\ x (t), which, as we shall show, is 
akin to a momentum. A self-consistent equation for this quantity can be found taking 
its derivative with respect to time: 

— x J xW = ~~ dt~^\ x ® + ^^U*) ( 23 ) 
In the previous expression, d t d x 9\ x ( t ) is the space derivative of f )20|) . Using fl22|) the two 
terms in d x 9\ x M cancel out, and fl23|) reads: 

= - d * D *Ut) +2cos(2d x d\ x(t) ) d x W x \ x(t) (24) 

This, together with eq. (122|) . provides a set of closed differential equations for x{t) and 
d x 9\ x ny The last thing we have to show is that the quantity d x 9\ x ^ is the average value 
of the momentum operator p, which reads: 

(p) = - / ^td x ip x 

x r (25) 

i.2 



= i \ \^ x \d x g(x,t) + / \r x \d x 9(x,t) 

J X J X 

Because of the form ( fT8l) of ip x (t) these integrals can be performed by the saddle point 
method (in the limit V — > oo). At the saddle point, d x g(x(t),t) = and d x 9(x(t),t) has 
a non zero value, thus (p) = d x 9\ x ( t ) +0(1/V). The same is true for (x) = x{t) + 0{l/V). 
Thus, rewriting (|2^|) and (I24p . one finds that the phase p(t) and the position x(t) obey 
the equations: 

^ = 4W X sin(p) = d p (D x - 2W X cos(2p)) 

A f \ (26) 
^ =-d x (D x -2W x cos(2p)) 

As a conclusion, the average values x(t) and p(t) of the position and momentum operator, 
which describe the global behavior of the wave packet, have the same time evolution 
as classical canonical variables with Hamiltonian H[x,p] = D x — 2W / x cos(2p). This 
property fits into the general picture of the propagation of wave packets in the semi- 
classical regime (see for example [31] for a mathematically oriented review) and, hence, 
was expected on general grounds as discussed previously. The generalization to more 
than one variable is straightforward. The width of the packet can also be evaluated. 
It is related to the separation in time of two classical trajectories initially separated by 
the width of the packet. When the effective dynamics is one dimensional, the classical 
Hamiltonian is integrable and periodic orbits separate linearly in time. In this case, 
since trajectories start from an initial distance 1/y/V, the typical time of separation 
is t ~ \/V. We check numerically that our analysis of the limit of large V is correct. 
In figure dj the exact quantum evolution, obtained by diagonalization of the discrete 
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5 10 15 20 25 

t 

Figure 1. The trajectory x(t) for U = 3.33, with initial conditions x(0) = 0.26 
of energy E = 0.2. The dots are obtained by numerical diagonalization of (|TU|) for 
V = 5000 with a sharp Gaussian initial condition. The line is the evolution according 
to the classical Hamiltonian. The two are equivalent on short times. Similar results 
hold for any trajectory and initial condition. 

equation ([7]), is compared to the classical evolution for short times. We find an excellent 
agreement. 

5.4- Different effective trajectories in the phase space 

We now study the effective classical dynamics by focusing on the phase space properties. 
Since the Hamiltonian f lT5|) is one-dimensional, and energy is conserved there are lots 
of constraints on the motion. In particular all the trajectories are integrable, and can 
only be periodic in x(t) except for separatrix trajectories. 

We find that the phase space is divided in two regions by a separatrix. To see this, one 
can characterize trajectories in terms of their turning points, and considering whether 
the momentum p is bounded or not. There are three types of turning points: the 
derivative % = AW (x) sin(2p) vanishes if either p = 0, p = tt/2 or W(2x') = 0. This 
last case is called "absorbing" for a reason that will be explained later. 

Let us compare three trajectories at different energies for U < U c , which are 
representative of all the cases encountered. For this purpose, two representations are 
shown: the evolution of x(t), pit) in figure [2], the phase space in figure [3] (right panel). 
The three types of trajectories are: 

• (A) has two p = turning points, and thus its momentum p is bounded. 

• (B) is a separatrix trajectory in the classical mechanics sense. The left turning 
point is at p = 0, and the right turning point at x = is "absorbing", the time 
taken to reach it (or escape from it) is infinite. In the equation of the motion, at 
the point W(x) = 0, the effective mass tend to infinity. 

• (C) has one turning point type at p = and one at p = n/2. p is growing infinitely 
large with time, which is not pathological because only p modulo ir enters the 
equations. 

The diagram in figure [3] (left panel) is similar in the spirit to the figures of effective 
potential shown to explain central motions in textbooks. It allows one to understand 
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Figure 3. a) Three trajectories A, B and C for U = 3.33 < U c . The energy of each 
trajectory, D + 2W and D — 2W as functions of x. b) Trajectories in the phase space, 
momentum 2p versus x. 



the dynamical evolution in a simple way. Actually, since the classical energy E = 
D(x) — 2W(x) cos(2p) is conserved, the values of x during the motion are restricted by 
the conditions D(x) - 2W(x) <E< D(x) + 2W(x) (note that W{x) ^ for x ^ 1, 0). 
Except for the separatrix, all turning points correspond to either p = 0&tE = D — 2W, 
or p = n/2 at E = D + 2W . Thus, all trajectories are delimited by the two turning 
points x a and Xf,, which are at the intersection between E and D(x) + 2W(x), and E 
and D(x) — 2W(x). The evolution of the system consists then in a periodic motion 
oscillating between x a and 

In the previous discussion we considered a specific form of the phase space and of 
D — 2W, D + 2W which are valid for certain values of U only. In the next section a 
careful investigation of the U dependence will be presented. However, before that, we 
want to stress that the effective potentials D — 2W, D + 2W can take three different 
qualitative form, depending on the coupling U. These are shown in figure HI The three 
regimes are separated by two special values of U: Ud and U c , see figure HI U c is the 
critical coupling of the quantum phase transition, and Ud = l/U c ~ 0.1715. The special 
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Figure 4. D(x) ± 2W(x), versus x for the three regimes of couplings U in the 
Bose-Hubbard model n™ ax < 2. Dark circles and squares are the superfiuid and 
Mott insulator ground states respectively. Dark dashed lines correspond to separatrix 
trajectories. 

regime U <Ud disappears when larger filling numbers n™ ax > 3 are included, so it is a 
peculiarity of the n™ ax = 2 case and, hence, not very relevant. 

6. Sudden quenches in the Bose-Hubbard model for n™ ax = 2 

We now discuss the dynamical evolution following a quantum quench and its dependence 
on the final and initial value of U. At t < the system is in the ground state at the 
coupling Ui. At t = the coupling is switched to Uf and the quench dynamics is 
computed for t > 0. In the following we call a quench "from superfiuid" when Ui < U c 
and "from Mott" when U > U c , and "to superfiuid" or "to Mott" when U f > U c and 
Uf < U c . The results of the following sections are summarized by the dynamical phase 
diagram shown in [6] (left panel) . 

6.1. Mott to superfiuid 

Starting from the Mott ground state x = 1, the trajectory is stuck at x(t) ~ 1 even 
for large times. In order to check this, let us linearize the equation of motion around 
x — 1. We use that x = dH/dp = 4W X sm(2p) and that sin(2p) can be extracted from 
E = D x — 2W x cos(2p) to obtain x = —^AW^ — ~D%. Thus, at dominant order in e we 
obtain the equation for e = 1 — x: 

e = e /r, r = 2/^(U c - U)(U - U d ) (27) 

The trajectory e(t) = is unstable, and since the wave function has a width 1/y/V, 
its typical evolution is given by e(t) = l/y/Ve 1 ^. Therefore, in the effective picture, 
the trajectory is stuck at the Mott ground state on times of the order of log(V). This 
result is a peculiarity, actually a pathology, of mean field models. Indeed, in a real finite 
dimensional system, spatial fluctuations drive the system away from the Mott state in 
a finite time. 
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6.2. Superfluid to Mott 

In a superfluid to Mott quench, the initial condition is given by the ground state packet 
characterized by {x = u ^ u ^ +1 ; p = o}. The trajectory after the quench is of the type (C), 
see quench Ql of figure [5j The superfluid order parameter |\Pq| (equation fl2])) oscillates. 

6.3. Superfluid to superfluid Ud < Uf < U c 

Depending on the value of C/j and Uf, there are three different types of dynamical 
evolutions in this case, that we call Q2, Q3 and Q4 in figure|5j Q2 is of type (C) and Q4 is 
of type (A). Q3, reached from a special Uf (dependent on Uf), corresponds to a singular 
separatrix of type (B) and of energy E = 0. In this case the packet relaxes exponentially 
to the Mott insulator ground state \x = 1), and concomitantly the superfluid order l^ol 2 
vanishes exponentially in time. This may be seen from the linearization (|27j), which 
becomes e = — e/r for a trajectory relaxing to x = 1. Approaching the transition 
oscillations take place on a time scale that diverges as — rlnQUf — Uf\). In consequence 
we find that a dynamical singularity, or transition, takes place at Uf. The values of Uf 
depends on U and thus defines a dynamical transition line Uf = (U c + Uf)/2, in the 
U,i,Uf plane. In figure [6] (right panel), as an example of singular behavior, we show the 
time average (|\l/o| 2 ) as a function of Uf for quenches starting from the non interacting 
case U = 0. 

It is interesting to compare to (|\I/o| 2 ) its equilibrium counterpart obtained from the 
microcanonical average corresponding to the same energy. From this we clearly see 
that the system is not thermalized. Actually, large quenches (Q2) have non vanishing 
oscillations of the superfluid order contrary to the equilibrium value which is instead 
zero (because it corresponds to an effective high temperature). Moreover, we find that 
the dynamical transition takes place at an energy at which the system, if it were relaxed, 
would be superfluid. Thus, the exponential relaxation of the superfluid order is a purely 
dynamical effect. 

6.4- Superfluid to superfluid Uf < Ud 

In this regime, there is a qualitatively new family of quenches Q6 (comparable to Q2). 
There is also a new special quench Q5 to the modified separatrix state at energy Ed, 
which gives rise to another dynamical transition. However, unlike the other quench 
cases, these effects are artifacts of the truncation n™ ax = 2 and disappear for n™ ax > 3. 
Therefore, they are not indicated in figure [6] (left panel), where the outcome of all 
possible sudden quenches is summarized. 
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Figure 5. Quench trajectories for the three regimes of Uf with D(x) ± 2W(x), like 
in figurelU Right panel (U c < Uf): In the quench Ql, at t = the packet state is at 
x = xq,p = (thus on the line D — 2W(Uf)), a position indicated by the arrow Ql. 
Then the trajectory of constant energy E is figured by the horizontal dark dash-dotted 
line. Left and center panel (Uf < Ud and Ud < Uf < U c ): Different initial conditions 
lead to qualitatively different quenches. The dark dashed lines are singular trajectories 
of infinite period. Quenches on these trajectories (Q3 and Q5) are at the dynamical 
transition. 




Figure 6. a) Dynamical phase diagram for the Bose- Hubbard model with n™ ax = 2. 
In the M area, within mean field, the system remains stuck to the Mott insulator 
ground state after the quench. Quenches from the supcrfluid phase are oscillating and 
similar to (A) or (C). The dynamical transition (B) separating the two is displayed 
as a dashed line, it meets the Mott phase at Uf = U c . b) Supcrfluid order (|*o| 2 ) 
as a function of Uf. Continuous line: time average after a quench. Dashed line: 
microcanonical average at the corresponding energy after the quench. The initial 
coupling is U\ = 0, but the evolution is qualitatively similar for all [/, with a dynamical 
transition. 



Dynamical transitions and quantum quenches in mean-field models 



18 



7. Bose-Hubbard model with weaker truncation n™ ax > 3 

7.1. Hamiltonian 

The previous analysis focused on two bosons or less per site. This constraint can be 
relaxed to include up to n™ ax number of bosons per site. When n™ ax is sufficiently high 
compared to the density, n™ ax 3> n, one recovers the behavior of the Bose-Hubbard 
model with no constraints on the occupation number per site. 

For a given maximum number of bosons per site, n™ ax , any symmetric wave function 
can be parametrized by the fractions Xi of sites with i bosons per site, i G [0, n™ ax \. Since 
the Xi are fractions, they verify Y2i x i = 1- Moreover the density n = ^2 i iXi is fixed, 
thus there are only n™ ax — 1 free variables left. The wave function is expanded in the 
symmetric basis like \ip) = ^2 x ipx(t)\x). The transition elements D(x) and W m {x) can 
be computed as done previously. For instance, there are 3 different types of transitions 
m when n™ ax = 3 and 6 when n™ ax = 4. Performing the classical equivalence for packet 
states, the resulting Hamiltonian can be put in the form 0: 



Specifically, when n™ ax = 3, if one chooses X\ and x 2 as free variables, the Hamiltonian 
is 

H(xi,x 2 ,pi,p 2 ) = D - 2WiCOs(pi +p 2 ) - 2W 2 cos(2pi - p 2 ) - 2W 3 cos{p 1 - 2p 2 ) 

W x = J(3x £ix 2 :r3) 1/2 W 2 = Jxi(2x x 2 ) 1/2 W 3 = Jx 2 (6xix 3 ) 1/2 

D = x 2 + 3x 3 - J(x xi + 2xix 2 + 3x 2 x 3 ) 

xo = 1 — xi — x 2 — X3 x 3 = |(n — xi — 2x 2 ) 

Notice that because < Xj < 1, there are constraints on the possible values of Xj, such 
as x 2 < (n — xi)/2. For all riJ£ ax there is a Mott insulator to superfluid quantum phase 
transition at some coupling U c if the density n is an integer (lower than n™ ax ). Above 
the critical coupling U > U c the ground state is a Mott insulator x n = l,x^ n = (all 
sites have n bosons), and below U < U c the ground state is superfluid with all Xj ^ 0. 

1.2. Regularity of trajectories after a quench 

In the previous case n™ ax = 2, the effective dynamics was one-dimensional, and thus 
integrable. For n™ ax > 2, the classical dynamics takes place in two or more dimensions, 
and the trajectories may be either regular or chaotic. In order to characterize them 
we study their regularity properties. For chaotic trajectories, neighboring trajectories 
separate exponentially in time in the phase space y = {xj,pj}, like 5y(t) ~ exp(Xt)5y(0). 
The rate of separation A is the largest Lyapunov exponent [32]. We computed A 
for n™ ax = 3 at density n — 1 for several trajectories, using a simplified version of 
the traditional Gram-Schmidt orthonormalization of the Lyapunov vectors. Roughly 
speaking, if we write the Hamiltonian evolution under the form y\ = fi(y), the deviation 
satisfies 5yi = djfi(y)8yj. This can be integrated numerically and normalized at each 
step to avoid an overflow (the orthogonalization step is dedicated to find all Lyapunov 
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Figure 7. a) Lyapunov exponents of the trajectories with initial conditions 
{xi,X2,Pi = p-2 = 0}, for U = 2.86 plotted in levels of gray. The bright zone is 
regular (numerically A < 0.05), the dark zones are chaotic. The Lyapunov exponents 
corresponding to initial conditions given by ground states obtained varying Ui are 
indicated by a continuous line. The dashed line indicates initial conditions with zero 
energy. The dynamical transition corresponds to their intersection, b) Dynamical 
phase diagram for n™ ax = 4. The dynamical transition is at the dashed line. The 
dotted line is the dynamical transition when n™ ax = 2, for comparison. 



exponents, whereas here we only need the largest one). The Lyapunov exponents for 
trajectories with different initial conditions {x\,x l 2 ,Pi = P2 = 0} are plotted in figure 
[7J in the superfluid phase U = 2.86 < U c . Trajectories are regular in some regions 
of the space (periodic or quasi periodic) with A = within the error bar, whereas 
some other regions are chaotic with A > 0.1. The quench from E = is exactly at 
the dynamical transition, the corresponding trajectory is chaotic. For U > U c , when 
the ground state is a Mott insulator, all trajectories are regular. We notice that the 
regularity of a trajectory affects the time of spreading of the packet (determined by the 
time of separation of two neighboring trajectories). Actually, the time of separation is 
typically polynomial t ~ V x l a for a regular motion but only t ~ In V for a chaotic one. 

7.3. Dynamical transition 

In quenches from the superfluid phase Ui < U c , like in the previous case r0£ ax = 2, a 
dynamical transition occurs at the special coupling Uj where the final energy equals 
the energy E — of the unstable Mott trajectory x n = l,x^ n = 0. Direct evidence 
of this transition is given by the singularity of the time averaged superfluid order l^ol 2 
as a function of Uf for a given C/j. In figure [8b, we compare this singularity in l^ol 2 
for n™ ax = {2,3,4,5} and find that the dependence in n™ ax of the divergence is very 
weak beyond n™ ax = 4. The two curves n™ ax = 4 and n™ ax = 5 are not distinguishable, 
because they are identical up to 0.01%. 

We show the dynamical phase diagram in figure [7b for n™ ax = 4 and unit filling factor. 
We observe that the transition line is near to the transition line for n™ ax = 2, and that 
they are asymptotically equal around U c . Because the probability of having more than 
4 bosons on the same site is extremely small, of the order of 0.01% when n™ ax S> 4, we 
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Figure 8. a) Evolution of x(t) and p(t) with time, n^ iax = 3, n = 1 and Ui = 1. Left 
panel, Uf = 2.5, and right panel, Uf = 3.29. The dynamical transition is at Uj = 3.21. 
The left panel is before (E < 0) the dynamical transition, the right one (E > 0) is 
after, b) Superfluid order parameter (l^ol 2 ) as a function of Uf for n = 1, £Zj = 3, 
with ?i™ aT = 2,3,4 and 5. The curve ri% lax = 2 is shifted of 0.025 along the U f axis 
for comparison. 



can safely assume that this phase diagram is quantitatively representative of the phase 
diagram without truncation. 

Even though the existence of the dynamical transition for any n™ ax is beyond 
doubt, because the singularity is numerically manifest, more precise results on this 
transition, even for n™ ax = 3, are hard to provide. Some features of the case n™ ax = 2 
persist ; for example, at the dynamical transition, the momentum 2pi — p 2 becomes 
unbounded, see Fig. [8^l. The fractions are also oscillating, but without definite 
period. They are either quasi-periodic for low A regions, or chaotic. Approaching 
the dynamical transition, the typical time of return to X\ ~ 1 is increasing, possibly 
logarithmically divergent in Uf — Uf as suggested by the numerical divergence in figure 
|8b. Unfortunately, the analysis of the singularity for n™ ax > 3 turns out to be out 
of reach. The numerical integration of classical equations of motions is exponentially 
sensible to numerical errors and thus not reliable. One could imagine that the trajectory 
at E = is an unstable manifold, which would support the existence of a singularity for 
trajectory arbitrarily close to it. However, it is hard to decide whether the trajectory 
at the dynamical transition goes arbitrarily close to the point X\ = 1 (ground state of 
the Mott insulator), at which the trajectory is exponentially slowed down. A possible 
scenario, in which the surface E = is ergodic, does not seem to be validated by 
numerical integration of trajectories. 
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8. Quenches in the generalized Jaynes-Cummings model for super-radiance 
transition 

8.1. Super-radiance quantum phase transition 

We now consider the generalized Jaynes-Cummings model (j4|) and derive its quench 
dynamics. For simplicity and without loss of generality we focus on Q = N/2. Thanks 
to the conservation of Q = S z + b'b, the spin degrees of freedom can be accounted for 
in terms of the bosonic ones. The states can be parametrized by the density of bosons 
n only, in the Fock basis |n) bosons <S> \S; m} 8pin : 

\n) = \Nn) hosons <g> \N/2; N(l/2 - n)) spin 

The action on \n) of some useful operators can be readily computed at the dominant 
order in TV 

S z \n) = N(l/2-n)\n) 
tib\n) = Nn\n) 

S + b\n) = V^N^/N 2 /A-m 2 \n - l/N) 
- n\n - l/N) 

From these equations, the transition elements D n = j^{n\H\n) and W n = ± 
\jN\H\n) are easy to compute, and yield an effective Hamiltonian on n and <fi its 
conjugate momentum H[n, (f>] = D n — 2W n cos(4>) 

H[n, (f)} = -2Xn + 2nVl - n cos(0) (30) 

which is the semi-classical approximation of the original Hamiltonian obtained in [25J, 
up to a shift in the phase <p — > (p + n. In [25], the authors studied sweeps from the empty 
state (n) = 0, and the link to the Landau-Zener problem. Here, we focus on sudden 
quenches and the related dynamical transition. The ground state of this effective system 
is at GS = Ti and 

{0 if A < — 1, Standard ground state 

| (3 — A 2 — a/ A 2 (3 + A 2 ) if — 1 < A < Super-radiant ground state 
| (3 — A 2 + a/ A 2 (3 + A 2 ) if A > Super-radiant ground state 

The super-radiance quantum phase transition is the result of a competition between 
the potential energy term and the kinetic energy term, which can be understood if the 
quantum Hamiltonian is written in the \n) basis: 

H = -2XNh + -^=NJn(l-h)(b ] + b) 
y/N V 

Qualitatively, the potential term — 2\b^b encourages filling for A > and discourages it 
for A < 0. The "kinetic" contribution ~ g(tf + b) can lower the energy. As a matter 
of fact, the sign of g is almost irrelevant. Changing the sign of g leads to W n — > —W n . 
However, the new Hamiltonian can be mapped into the old one by the translation 
p — > p + 7r. In consequence, if the ground state of a given Hamiltonian is at p = 0, 
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Figure 9. a) Ground state n GS as a function of the parameter A. The quantum phase 
transition takes place at A = — 1. b) Quench diagram for all A.; — > Xf. The quench 
trajectory is of bounded momentum <f> (light gray) or unbounded (white). On the 
dashed line, the separatrix states are not singular (absorbing). The thick line is the 
dynamical transition, where the trajectories are singular. For A^ < —1, in the dark gray 
region, trajectories are definitely stuck at n = 0. Otherwise, there is a slow relaxation 
on times of order log N. 



ip n ~ exp[— (n — n GS ) / (a / \/N) 2 ], then the ground state of the Hamiltonian with the 
opposite sign of g is at p = tt, xjj n ~ exp[— iNnn — (n — n GS ) / (cr / y/N) 2 ], identical but 
with a staggered sign. 

8.2. Sudden quenches 

Previously, we have shown that the characterization of turning points plays an important 
role for understanding the dynamical behavior. In the generalized Jaynes-Cummings 
model, there are two special turning points, n c = {0, 1}, identified by the equation 
^jf = W n sm((j)) = 0. Like in the Bose-Hubbard model, only one gives rise to singular 
trajectories. This is readily seen, for example, computing the time taken to reach the 
singular point and having zero kinetic energy at this point: 

T= dn ^ = / dn(2W n sin(0))- 1 

Using the conservation of the energy E = D n — 2W n cos(0), we can eliminate sin(0) = 
[1 - (E(n c ) - D n ) 2 /(4iy n 2 )] 1 /2 an d obtain 

/n c 
dn[AW 2 n - (E(n c ) - D^y 1 ' 2 

Around n = 1 — e, the time is finite T ~ j Q de l/v^e, whereas around n = e, the time 
T ~ Jo de l/v^ diverges. In other words, the trajectories touching the turning point 
n = 1 have a finite period, whereas the trajectories touching n = have infinite period. 
The situation is anologous to the one already studied for the BHM: the absorbing state 
n c = plays the same role of the Mott state, and gives rise to the dynamical transition. 
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To identify the critical value Aj at which the transition takes place, one can compute 
the energy after a quench A j — > Xf and compare it to the energy of the system at rest 
in n = 0. One can check that a trajectory n(t) starting with the same energy, actually 
energy zero, is indeed driven to the absorbing point n c = 0. This dynamical transition 
occurs on the line A^ = — a/1 — n GS (Xi). A linearization similar to (1271) shows that the 
relaxation is exponential at the transition with a characteristic time r" 1 = 2a/1 — A 2 . 
Around the transition, the period of a trajectory diverges as r\n(\Xf — Aj|). 

The dynamical phase diagram is shown in figure [9b. For Aj > A c , the initial state 
is in the broken symmetry phase {n(t = 0)) ^ 0. The situation is very similar to the 
BHM one: the filling number n(t) oscillates in time after the quench, and there is one 
region with bounded effective momentum <p, and two regions with unbounded <fi. Notice 
that there are two separatrices between the three regions, and that one is not singular, 
and does not give rise to a dynamical transition, whereas the other is. The singular one 
is the one for which the absorbing state n = is met, corresponding to the restored 
symmetry state. For quenches with Aj < A c , the initial state is empty (n(t = 0)) = 0, 
and either drifts away from zero on large times (of the order of logiV) if Xf < 1, and 
otherwise, the state is definitely stuck at n — 0. The bounded or unbounded nature of 
the final state is also indicated. 

Finally, we remark that eigenstates of the completely connected model can be 



written within a WKB expansion, this is done for completeness in Appendix B 



9. Quenches in the Ising model in a transverse field 

In our last example, let us consider the dynamics due to quantum quenches in the 
transverse field Ising model. To derive the classical effective dynamics, one can use the 
fact that the large spin limit is also the classical limit, see e.g. [21]. To do so, one 

should make the substitution S — > S = S{sin 9 cos <p, sin 9 sin <p, cos 9}. The kinetic term 
(needed to describe classical dynamics) in the classical Hamiltonian of the large spin 
limit can be derived using path integral for spins [33]. Here we instead use our generic 
method based on site permutation symmetry. We denote spin symmetric states of total 
momentum S = N/2 with the notation \s). They are such that S z \s) = Ns\s) with 

s G [-1/2, 1/2]. At dominant order in N, S x \s) = N^J\ - s 2 (\s + ±) + \ 8 - £)) /2. 
One can proceed as usual to get the effective Hamiltonian: 

D s = —(s\H\s) = --s 2 

N 2 



i i - r i 

W s = (s+ — Hs) =-J--s 2 

iV V N l 1 1 2 V 4 

H[s,(f)} = D s - 2W S cos = --s 2 -rJ--s 2 cos( 



As expected, the Hamiltonian is very reminiscent of the classical Hamiltonian for a 
single rotor once the change of variables s — > cos 9 is made. Indeed, the equations of 
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Figure 10. D(s) ± 2W(s) for r = 0.3 < T c and r = 0.6 > T c . The dark square is 
the quantum paramagnetic ground state, the dark circles are the ferromagnetic ground 
states. The dark dashed line is the separatrix, quenches on this trajectory are at the 
dynamical transition. 



motion obtained here are the same as is [21]. In this paper, the authors studied the 
quench from the paramagnetic phase to the ferromagnetic phase and AC dynamics. 

Proceeding as for the BHM, we analyze the phase space of symmetric trajectories. 
This is different in the ferromagnetic and paramagnetic phases, see figure [TO) In the 
latter, obtained for T < T c , there is a separatrix. For sudden quenches from the 
ferromagnetic phase, where s ^ 0, to final values of the transverse field Tf < T c , there is 
a dynamical transition at Tj = | + Tf. At the transition, the trajectory is a separatrix, 
and s decreases exponentially in time to s — 0, which is the quantum paramagnetic 
ground state. The symmetry is dynamically restored at this point but also for larger 
quenches: contrary to the two previous examples, for quenches beyond the dynamical 
transition Tf > Vj, the trajectories are symmetric in s — > —s, and the magnetization is 

oscillating around zero. Averaging over time, the order parameter is zero (S z ) = 0. 
10. Discussion 

Before concluding, two points need to be addressed further: (1) the possible link between 
dynamical transitions and equilibrium quantum phase transitions and (2) the pro and 
cons of our approach concerning the physics of finite dimensional systems. 

• In three examples, we found that dynamical transitions occur in systems 
characterized by a quantum equilibrium phase transition. Therefore, despite the 
fact that the dynamical transition is not directly related to the equilibrium one, as 
we stressed before, it is natural to wonder whether instead its existence is related 
to it. In support of this, if one studies the completely connected Bose-Hubbard 
model in a parameter regime (non unitary fillings) where there is no equilibrium 
quantum phase transition then one finds that there is no dynamical transition 
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either. Actually, we found the very same phenomenon in all models we studied, 
and this was also noticed for the Fermionic Hubbard model [TT]. We certainly 
cannot address this issue in general, also because it is still far from clear what "the 
dynamical transition" is beyond mean-field theory. We shall instead endeavor to 
provide a partial answer at the mean-field level. Generically, we observe that in 
all models we studied, there is a regime in which the symmetric state (the Mott 
insulator, the paramagnet, etc.) is not the ground state but it plays the role of an 
excited metastable state that can trap the system for an infinite time. This gives 
rise to the dynamical transition. If the effective motion is one dimensional, this state 
is absorbing and belongs to a separatrix. If the effective motion has two or more 
degrees of freedom, this state seems to belong to an unstable manifold. Thus, it 
appears that the singularity at the dynamical transition is a result of this singular 
trajectory, which in turn is a consequence of the existence of a quantum phase 
transition. A physical argument to support this idea is the following. For a quench 
from the unbroken symmetry phase to the broken symmetry phase in completely 
connected model, the dynamical evolution is expected to be very slow. The reason 
is that the symmetry must break up and the order must develop starting from this 
point. This, however, is a local phenomenon, for example due to domain growth in 
finite dimensions, and therefore generically absent in completely connected models, 
except on times diverging with the system size: indeed, in our examples, we find 
that the typical time scale to depart from the initial value x(t = 0) is of the 
order of log V. The point is that this slow trajectory is precisely the time reversal 
of the singular trajectory that is responsible for the dynamical transition. To 
summarize, the existence of a quantum phase transition implies the existence of 
slow trajectories. They correspond to the escape from the symmetry unbroken state. 
Their time reversed counterparts are the trajectories responsible for the dynamical 
transition. Thus, the existence of a phase transition indeed appears to be 
intertwined with the existence of a dynamical transition within mean- 
field. 

Let us now discuss the range of applicability of our results for finite-dimensional 
systems. First, our analysis can be put on the same footing as the Gutzwiller 
Ansatz and functional integral saddle point approximations, since the corresponding 
dynamical equations of motion coincide (see Appendix A for the former equivalence 
and also [22] ). Clearly, our mean-field approximation misses several essential 
physical effects. Relaxation and thermalization, for which spatial and temporal 
fluctuations must be taken into account, are not present. Furthermore, the 
dynamics of a quench from the unbroken symmetry phase to the broken symmetry 
phase is not properly described for the reasons mentioned above. Inhomogeneities, 
topological defects, domain growth are out of reach. 

On the contrary, our mean field approximation is expected to capture well the 
evolution of local quantities at least on short times, as long as they are homogeneous 
across the system. Contrary to many other approaches, it is not perturbative in 



Dynamical transitions and quantum quenches in mean-field models 



26 



parameters (U, J, T . . . ) of the Hamiltonian. Unlike pertubative approaches, 
it can describe the existence of dynamical transitions, and phenomena 
related to global observables. 

11. Summary and outlook 

The general purpose of this article was twofold. First, we described a method to 
study the quantum quench dynamics of generic completely connected models, providing 
a mapping to effective classical dynamics. There are two reasons for considering 
completely connected models. In some cases, such as in the generalized Jaynes- 
Cummings model, and the Dicke model, they provide the correct physical description. In 
others, such as the transverse field Ising model or the Bose-Hubbard model, they lead to 
an approximative treatment of finite dimensional systems. In the latter case, the range 
of validity of the approximation is possibly limited to short times only. Our second aim 
was to study and to reveal the existence of out of equilibrium dynamical transitions 
induced by quantum quenches. In agreement with other works [HI [13] we showed that 
within the mean-field approximation dynamical transitions occur in quenches from the 
broken symmetry phase to other regions of the broken symmetry phase and that this 
is a quite generic phenomenon for systems that display a quantum phase transition at 
equilibrium. 

Clearly, a main and pressing question is to understand what this transition is really, 
beyond mean-field theory. Does it become a cross-over in finite dimensions? If it remains 
a bona fide transition, what are its critical properties? In the d = 1 BHM, a non- 
monotonic behavior of the propagation velocity as a function of Uf has been found in 
[34] by t-DMRG. This could well be a signature of a cross-over related to the dynamical 
transition found in mean-field. An exact diagonalization study of hard core bosons in 
d = 2 suggests the existence of a cross-over or a singularity of the revival time after 
a quantum quench [23]. From the analytic point of view going beyond mean- field is a 
difficult task. An argument for the occurrence of dynamical transitions beyond mean 
field is given in [13], within an imaginary time path integral formalism. Clearly, it is 
crucial to take into account fluctuations not captured by mean-field theory. This was 
started to be done in [35J. This work indeed suggests that these fluctuations could 
alter substantially the mean-field results. Promising ways to capture fluctuations are 
the projector operator formalism [36], 1/z expansions [37] (z is the connectivity of the 
lattice). Field-theoretic methods are also being developed, as in [19] for the Bose- 
Hubbard model in the grand canonical ensemble; the use of two particle irreducible 
actions also seems promising [38] . 
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Appendix A. Gutzwiller Ansatz 

The previous analysis, based on the symmetry between sites, seems to be related to 
completely connected models only. It is instead much more robust than it looks, since 
it is actually equivalent to the time dependent Gutzwiller Ansatz. The latter can be 
formulated as a mean field approximation for a finite dimensional model: it consists in 
providing a trial wave function dependent on some physical parameters that are chosen 
through a maximization procedure. In its equilibrium version, the estimated ground 
state is simply found minimizing the energy. The dynamical formulation is a bit more 
elaborated. For clarity, we show here how to build it for the case n™ ax = 2 only. The 
following is an adapted version of the formulation developed in [TT] . 

Let us define first the "full state" \F) = |0) + |1) + |2), and let Pj be the projector 
on site / on the state \j), where j G {0, 1, 2}. The Gutzwiller Ansatz wave function is a 
trial wave function with 6 real, time dependent parameters Xj(t), Pj(t): 

I^ga) = n e'toWA-H**,) (v^P + V*TA + V^A) \F) 

(Al) 



f] (e^v^|0) + e^v^Ill) + e*>^\2)) 



Note that for simplicity we dropped the site index. The x n are the projector average: 
x n = {^}\P n \ij}) . As we shall see, they will turn out to coincide with the fraction of sites 
with n bosons. Obviously, because this is only a trial wave function, it can not satisfy 
exactly the time evolution: i.e. ih^\ip GA ) ^ H\ip GA ). However, one can enforce it in an 
approximate way. This allows one to obtain the time evolution of Xj(t) and Pj(t). We 
impose the following constraints: 

• The projectors in the Heisenberg representation satisfy on average the 
Heisenberg equation of evolution Pf 1 = i[H,P^]. This amounts to computing 
the average of P in the Schrodinger picture xi(t) = (i[)\Pi\ij;) using 

Xl =MAM + MA|^> (A2) 

Pi,HV- 



• The energy is conserved E = (ip\H\x/j). 

Let us call W[xi,pi\ = (if)\H\i/)) the energy as a function of the parameters. The evolution 
(1A.2I) can also be written, using the explicit form of \ip) (lA.ljl : 

x x =-ilM?\Mi>) = ^^>\KW) 

= &H[x ltPl ] Pl ^ 

dpi 

This evolution is very reminiscent of the first Hamilton equation of motion. Furthermore, 
the conservation of energy E = dn ^ Pl ^ pi + —j^— xi = leads to the second Hamilton 
equation: 

= _mp I i 

dxi 
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As a consequence, the Gutzwiller Ansatz wave function is indeed equivalent to 
an effective classical Hamiltonian evolution of the variables xi and pi. The effective 
Hamiltonian H[xi,pi] = (ip\H\ip)[xupi\ must be computed for the model at hand, and is 
indeed the same as the one discussed in section HI For instance, for the Bose-Hubbard 
model with truncation n™ ax = 2, one recovers (115j) . 



Appendix B. Eigenstates in completely connected models 

Using a Wentzel-Kramers-Brillouin (WKB) approximation [39], eigenstates of the 
Schrodinger equation ( {TBI) can be found in the limit of large V. To do so, we write 
the eigenstates 4>e{%) within the WKB approximation. The derivation is standard, one 
looks for stationary solutions of the form ip(x) = A(x)e tVS ^ of the equation (113j) . at 
dominant order in V . One finds 

<p E {x) = C exp (±iV [ dx'p(x')] (B.l) 



2W(x)^sm(2p(x)) 

where C is a normalization constant, E is the energy, p(x) satisfies the implicit equation 
E = D x — 2W X cos(2p(x)), and x is the left turning point of the classical trajectory of 
energy E. 

Using this expression, we can draw a parallel between the eigenstate and the 

effective classical trajectories of energy E. For an observable f(x), the average reads 

(faimife) =\c\ 2 [ dx 



J x 4PU(x) sm{2p{x)) , B 2 s 

= fat L dx t m =W) 

In the second line, we refer to the effective classical motion ^| = and the average 
f(x) is the average over the effective classical trajectory of f(x) over one peridod. In 
words, we found that in the semi classical regime, the average over a quantum stationary 
state of energy E is given by the average over one period of the classical trajectory of 
same energy E. Although this result is not often mentioned, this is a direct consequence 
of the WKB expression of the wave function. 
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